\defmodule{KorobovLatticeSequence}

\pierre{This class is not yet fully implemented}

This class implements Korobov lattice sequences, defined as follows.
One selects a \emph{basis} $b$ and a (large) multiplier $a$.
For each integer $k\ge 0$, we may consider the
$n$-point Korobov lattice with modulus $n = b^k$ and multiplier
$\mbox{\~a} = a \bmod n$.
Its points have the form
\eq
 \mathbf{u}_i = (a^i (1, a, a^2, \ldots) \bmod n) / n
       = (\mbox{\~a}^i (1, \mbox{\~a}, \mbox{\~a}^2, \ldots) \bmod n) / n
                                        \eqlabel{eq:Korobov-seq1}
\endeq
for $i=0,\dots,n-1$.
For $k = 0,1,\dots$, we have an increasing sequence of lattices
contained in one another.

These embedded lattices contain an infinite sequence of points that 
can be enumerated as follows \cite{vHIC01a}:
\eq
 \mathbf{u}_i = \psi_b(i) \left(1, a, a^2, \ldots \right) \bmod 1.
                                        \eqlabel{eq:Korobov-seq2}
\endeq
where $\psi_b(i)$ is the radical inverse function in base $b$,
defined in \class{RadicalInverse}.
The first $n=b^k$ points in this sequence are exactly the same as
the $n$ points in (\ref{eq:Korobov-seq1}), for each $k\ge 0$.

\bigskip\hrule\bigskip
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{code}

\begin{hide}
/*
 * Class:        KorobovLatticeSequence
 * Description:  
 * Environment:  Java
 * Software:     SSJ 
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author       
 * @since

 * SSJ is free software: you can redistribute it and/or modify it under
 * the terms of the GNU General Public License (GPL) as published by the
 * Free Software Foundation, either version 3 of the License, or
 * any later version.

 * SSJ is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.

 * A copy of the GNU General Public License is available at
   <a href="http://www.gnu.org/licenses">GPL licence site</a>.
 */
\end{hide}
package umontreal.iro.lecuyer.hups;


public class KorobovLatticeSequence extends KorobovLattice \begin{hide} { 
   int base;         // Base for radical inversion
   int inverse;      // global variables for radical inverssion,
   int n;            // since bloody JAVA cannot pass references

   // Method modPower is inherited from Rank1Lattice.
\end{hide}
\end{code}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection*{Constructor}

\begin{code}

   public KorobovLatticeSequence (int b, int a) \begin{hide} {
// Pas termine: ne fonctionne pas
      super (2, 3, 1);
      if (a < 1)
         throw new IllegalArgumentException
             ("KorobovLatticeSequence:   Multiplier a must be >= 1");
//      dim       = Integer.MAX_VALUE;
//      numPoints = Integer.MAX_VALUE;
      base = b;
throw new UnsupportedOperationException ("NOT FINISHED");
   } \end{hide}
\end{code}
 \begin{tabb}
  Constructs a new lattice sequence with base \texttt{b} and 
 \texttt{generator} $ = a$.
 \end{tabb}
\begin{htmlonly}
   \param{b}{number of points (modulus) is a power of b}
   \param{a}{multiplier $a$ of this lattice sequence}
\end{htmlonly}
\begin{code}\begin{hide} 

   // A very inefficient way of generating the points!
   public double getCoordinate (int i, int j) {
      int n;
      int inverse;
      if (i == 0)
         return 0.0;
      else if (j == 0)
         return radicalInverse (base, i);
      else {
         // integerRadicalInverse (i);
         n = 1;
         for (inverse = 0; i > 0; i /= base) {
            inverse = inverse * base + (i % base);
            n *= base;
         }
         return (double) ((inverse * modPower (genA, j, n)) % n) / (double) n;
      }
   }

   // ... has been unrolled in getCoordinate.
   private void integerRadicalInverse (int i) {
      // Attention: returns results in variables n and inverse.
      n = 1;
      for (inverse = 0; i > 0; i /= base) {
         inverse = inverse * base + (i % base);
         n *= base;
      }
   }
 
}\end{hide}
\end{code}
